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Abstract 

MicroRNAs (miRNAs) are small RNA molecules involved in the regulation of mammalian gene expression. Together with other 
transcription regulators, miRNAs modulate the expression of genes and thereby potentially contribute to tissue and species 
diversity. To identify miRNAs that are differentially expressed between tissues and/or species, and the genes regulated by 
these, we have quantified expression of miRNAs and messenger RNAs in five tissues from multiple human, chimpanzee, and 
rhesus macaque individuals using high-throughput sequencing. The breadth of this tissue and species data allows us to show 
that downregulation of target genes by miRNAs is more pronounced between tissues than between species and that 
downregulation is more pronounced for genes with fewer binding sites for expressed miRNAs. Intriguingly, we find that 
tissue- and species-specific miRNAs target transcription factor genes (TFs) significantly more often than expected. Through 
their regulatory effect on transcription factors, miRNAs may therefore exert an indirect influence on a larger proportion of 
genes than previously thought. 
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Introduction 

MicroRNAs (miRNAs) are small RNA molecules (-18 to 24 
nt) that regulate gene expression posttranscriptionally via 
messenger RNA (mRNA) destabilization or translational 
repression (Lee et al. 1993; Wightman et al. 1993; Lim 
et al. 2005; Friedman et al. 2009). MiRNAs are characterized 
by high sequence conservation across highly divergent spe- 
cies (Pasquinelli et al. 2000; Hertel et al. 2006; Chen and 
Rajewsky 2007; Stark et al. 2007), which reflects the effects 
of purifying selection due to their evolutionary importance 
as regulatory molecules. Expression profiling of miRNAs has 
been performed using cloning and Sanger sequencing 
(Landgraf et al. 2007), microarrays (Lim et al. 2005), and, 
most recently, by direct high-throughput sequencing 
(Creighton et al. 2009). High-throughput sequencing 
technologies have facilitated both the identification and 
the large-scale expression profiling of miRNAs in different 
tissues and cells. MiRNA profiling was also used to study 



expression across various developmental stages (Somel 
et al. 2010, 201 1) and in comparison between healthy 
and diseased tissue (Erson and Petty 2008). 

The destabilization of target mRNAs is the predominant 
mechanism of expression regulation by miRNAs (Guo et al. 
2010). It is therefore possible to assess the extent to which 
miRNAs shape gene expression by quantifying both miRNA 
and mRNA expression in the same samples. The correlation 
between the expression of single miRNAs and their target 
genes has previously been studied in cell culture by knocking 
out/down endogenously expressed miRNAs or by introduc- 
ing miRNAs into a specific target cell (Baek et al. 2008; 
Selbach et al. 2008; Guo et al. 2010). Genome-wide 
patterns of correlation between expression of miRNAs 
and their target genes in multiple tissues have been less 
widely explored. 

In primates, the main focus has been on discovery and 
annotation of miRNAs (Berezikov et al. 2005, 2006a, 
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2006b; Kawaji et al. 2008; Yue et al. 2008; Baev et al. 2009; 
Brameier 2010), but much less is known about expression 
variation within and between species, and how these pat- 
terns vary in different cell types and tissues. Three studies 
have correlated the expression of miRNA, mRNA, and pro- 
teins in prefrontal cortex comparing humans and rhesus 
macaques (Somel et al. 2010) and between humans, chim- 
panzees, and rhesus macaques (Somel et al. 201 1 ; Hu et al. 
201 1). These studies showed that miRNAs play a role in pri- 
mate development and aging (Somel et al. 201 0, 201 1 ) and 
that miRNAs with human-specific expression patterns target 
genes involved in neuronal functions (Hu et al. 201 1). How- 
ever, an understanding of the evolution of gene expression 
regulation in primates across multiple tissues has not yet 
emerged. This is of particular interest given that it has been 
hypothesized that the phenotypic differences between spe- 
cies may be better explained by changes in gene expression 
than by changes in DNA sequence (King and Wilson 1 975). 
Studies exploring differences in gene expression in multiple 
primate tissues have yielded lists of genes, which are differ- 
entially expressed in closely related species (Khaitovich et al. 
2005), but it is less clear which regulatory changes drive 
these differences. 

We have characterized miRNAs from multiple individuals 
in five different tissues (brain, heart, kidney, liver, and testis) 
from three primate species (humans, chimpanzees, and rhe- 
sus macaques) by lllumina sequencing. We examined the ef- 
fects of miRNA-mediated regulation by comparing with 
a gene expression data set generated from the same sam- 
ples as used for the miRNA expression profiling. Using these 
data sets, we address the following questions: 1) How big 
are miRNA expression differences between species and tis- 
sues? 2) Is there a particular group of genes regulated by 
miRNAs, which are differentially expressed between species 
and/or tissues? 3) Can the functional relationship between 
miRNAs and their target genes be determined by measuring 
their expression simultaneously? 



Materials and Methods 

Samples Library Preparation, Sequencing, and Base 
Calling 

All the individuals used in this study were adult males and 
suffered sudden death that did not involve the tissues sam- 
pled. A description of the samples is available in supplemen- 
tary tables S1 and S2 (Supplementary Material online). 

Total RNA was prepared as described in the lllumina Inc. 
manual "Small RNA Sample Preparation Guide" (Part # 
1004239 Rev. A lllumina Inc., San Diego, CA). lllumina Ge- 
nome Analyzer I and II sequencing runs were analyzed starting 
from raw intensities. A detailed summary about the platform 
a sample was sequenced on, how many cycles and 
which chemistry were used can be found in supplementary 



table S2 (Supplementary Material online). Base calling and 
quality score calculation were performed for all runs using 
the Improved Base Identification System base caller (Kircher 
et al. 2009), trained on <pX174 control reads of a dedicated 
control lane in each run. Reads with sequence entropy higher 
than 0.3 (sequence complexity filtering) and without any bases 
below a quality score of 10 were kept for downstream analysis. 

Composition of Samples and Mapping of Small RNA 
Reads 

We quantified expression of previously annotated miRNAs 
from miRBase (Griffiths-Jones et al. 2008; www.mirbase. 
org, release 15) for human, chimpanzee, and rhesus ma- 
caque. We mapped the sequenced reads to the official miRNA 
repository and the corresponding species' annotated 
miRNAs using PatMaN (Prufer et al. 2008) allowing zero mis- 
matches. Only reads with a length greater than 1 8 nt were 
used. The mature sequences were used as reference se- 
quences for each miRNA. If a read was a substring of 
a miRNA sequence or vice versa, this read was assigned 
to this miRNA. If a read mapped to multiple miRNA sequen- 
ces, the counts of this read were equally distributed to all 
matched miRNAs. 

To classify reads that did not map to miRNAs in miRBase, 
we mapped the remaining reads to multiple databases in 
order to distinguish alternative read sources. We used the 
gene database provided by Biomart (ensembl.org/biomart) 
using all ENSEMBL genes (version 59) with their untranscrip- 
ted sequence for human. To identify different structural 
RNAs, we used sequence annotation from UCSC of RNA 
genes based on the NCBI 36.1 (hg 1 8) human reference ge- 
nome. We obtained the sequences of piwiRNAs (piRNAs) for 
human from RNAdb (http://jsm-research.imb.uq.edu.au/ 
rnadb/). Sequences not mapping to any of these databases 
were aligned to the corresponding species genomes (hg 1 9, 
panTro2, rheMac2). The reads were distributed to the differ- 
ent categories in the following hierarchical order: miRNA, 
piRNA, rnaGenes, genes, unknown but mapped to the 
genome and not mapped to the genome. 

mRNA Expression Data 

We used mRNA sequence tags quantified using the lllumina 
NG III Digital Gene Expression approach (Velculescu et al. 
1995; Kircher 2011). Samples were obtained from brain 
(prefrontal cortex), heart, kidney, liver, and testis tissues 
of male humans, chimpanzees, and rhesus macaques. 
The samples overlapped extensively with the samples used 
to generate the miRNA expression set (supplementary 
table S1, Supplementary Material online). 

Normalization of miRNA and mRNA Data 

We normalized both the miRNA and the mRNA data using 
the R package DESeq (Anders and Huber 2010). The 
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package performs a negative binomial fit to the samples. In 
addition, we log transformed the normalized data. 

Expressed miRNAs/Genes and Expression Differences 

We defined a miRNA as expressed if more than one tran- 
script was sequenced in one species and one tissue. A gene 
was considered to be expressed as defined in Kircher et al. 
(201 1 ), that is, all used samples had to have of at least three 
mapped reads. Expression differences were calculated using 
the R library DESeq. The function nbinomTest provides a P 
value for each compared miRNA. In addition, a correction 
for multiple testing based on the Benjamini-Hochberg pro- 
cedure was performed. We defined all miRNAs with an ad- 
justed P value < 0.05 as differentially expressed. The same 
procedure was performed for the mRNA data. 

Expression of Non-annotated miRNAs 

A miRNA was defined to be expressed in another species 
than human if in at least one tissue the expression level 
was bigger than zero for at least four individuals and 
a Wilcoxon rank sum test showed no significant (alpha = 
0.05) difference between the human samples expression 
values and the respective species samples expression values. 

Species and Tissue Effect 

To determine the effect of species and tissues in our data set, 
we used the set of miRNAs with at least ten transcripts in 
each species for the analyzed tissue. Additionally, only 
miRNAs annotated in all three species were included. We 
used the normalized data set. The fraction of variance 
explained by tissues and species was calculated by comput- 
ing the fractions of sum of squares explained by the factor 
tissue and the factor species in a linear model. 

miRNA Sequence Evolution 

We calculated miRNA average conservation based on scores 
obtained from the multiple alignments of 46 vertebrate spe- 
cies using phastCons46way (Siepel etal. 2005). For humans, 
we calculated miRNA single nucleotide polymorphism (SNP) 
density using all SNPs from dbSNP (v. 135) (Sherry et al. 
2001). 

Newly Predicted miRNAs 

We used miRDeep (Friedlander et al. 2008) to identify po- 
tential new miRNAs. For this purpose, we merged reads of 
all samples in a species. Predictions with positive miRDeep 
scores and in orthologous regions (UCSC liftOver; Hinrichs 
et al. 2006) of all species were used for further investiga- 
tions. Newly predicted miRNAs that were found in ortholo- 
gous genomic regions of all three species were submitted to 
miRBase. Accession numbers form miRBase are assigned 
after publication acceptance. 



Target Prediction 

We obtained miRNA-binding sites for all mRNA genes from 
the TargetScanS (Lewis et al. 2003) database. The predictions 
included both conserved- and nonconserved-binding sites. 

Functional Gene Ontology Analysis 

We used the Gene Ontology (GO) (Ashburner et al. 2000) 
and the hypergeometric test from FUNC (Prufer et al. 2007) 
to test for enriched GO categories among gene groups. 
Genes that were regulated by at least one miRNA with dif- 
ferent expression were coded "1," whereas all genes tar- 
geted by miRNAs that were expressed and showed no 
differential expression got the label "0." FUNC-hyper was 
run with parameter -c 10, which includes categories with 
at least ten genes in it. As our significance measure (alpha = 
0.05), we used the familywise error rate. 

Random Distribution Calculation 

Significance levels were computed by calculating a random 
distribution and comparing the observed values to this dis- 
tribution. The random distribution is computed by randomly 
assigning expressed miRNAs to expressed genes from the 
prediction list. Each miRNA and gene had the same chance 
of assignment. For the correlation between tissues depen- 
dent on the number of binding sites for expressed miRNAs, 
we randomly permuted the number of binding sites be- 
tween the miRNA and the mRNA pairs. The random assign- 
ments were performed 1,000 times. 

Results 

Small RNA Composition and Annotation 

We sequenced 73 small RNA libraries (25 from humans, 24 
from chimpanzees, and 24 from rhesus macaques) derived 
from 5 individuals of each of the species (human, chimpan- 
zee, and rhesus macaque) for the tissues brain, liver, and 
heart. Five individuals of each species were sequenced for 
kidney (except chimpanzee with four individuals) and testis 
(with four rhesus macaques). When we mapped the reads to 
miRBase (www.mirbase.org, release 15; Griffiths-Jones 
et al. 2008; supplementary table S1, Supplementary Mate- 
rial online), the majority of reads matched known miRNAs 
(median 59%). However, testis yielded a consistently smaller 
fraction of matches for all species (median 9%). A consider- 
able proportion of reads in testis were assigned to piRNAs, 
which are a different class of small RNAs that are mainly 
found expressed in the germ line (Girard et al. 2006; Aravin 
et al. 2007) (fig. 1a and b). An excess of 5' uridine (U) res- 
idues is a signature of the piRNA family (Malone and Hannon 
2009). This 5' uridine (U) enrichment in piRNAs is clearly vis- 
ible in testis-derived molecules (55% U in first position) com- 
pared with molecules from other tissues (supplementary 
fig. S1, Supplementary Material online). In addition, we 
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Fig. 1. — Small RNA libraries composition and variance distribution for miRNA expression, (a and b) Small RNA composition for a sample of human 
brain and human testis, respectively, (c) Percentage of miRNA expression variance explained by factors tissue and species and their interaction, 
(d) Percentage of miRNA expression variance explained by different factors in tissues. 



found that in testis a higher fraction of reads failed to be 
assigned to small RNAs, structural RNAs, or mRNAs. How- 
ever, these reads could be mapped to the respective genome 
and showed an even higher fraction of uridine in the first 
position than observed for the reads mapping to annotated 
piRNAs suggesting that they represent unannotated piRNAs 
(supplementary fig. S1, Supplementary Material online). 

Combining the data of all five tissues, we identified 585 
of the 718 (81%) known miRNAs in human, 431 of the 530 
(81 %) known miRNAs in chimpanzee, and 399 of the 502 
(79%) known miRNAs in rhesus macaque (table 1). Al- 
though the number of different miRNAs detected is similar 
between tissues, brain and testis showed consistently more 
expressed miRNAs than other tissues (table 1). 

When comparing the repertoire of miRNAs between spe- 
cies, we found some instances in which a particular miRNA 
was annotated in human but not in one or both of the other 
species. Using our rhesus and chimpanzee expression data, 
we were able to detect 1 6 of these unannotated miRNAs in 
chimpanzee and 27 miRNAs in rhesus macaque. The expres- 
sion level of these miRNAs in these other species was gen- 
erally comparable with the expression level in human. 

In order to identify miRNAs independently of their 
miRBase annotation, we applied the miRDeep algorithm 
(Friedlander et al. 2008) to our data set. miRDeep identifies 



miRNAs from deep sequencing of small RNA libraries by 
comparing the position and frequency of reads with the sec- 
ondary structure of the predicted pre-miRNA. The algorithm 
provides a combined score indicating the reliability of the 
prediction; the more positive the score the more reliable 
the prediction. To reduce the false positive prediction rate, 
we applied a score cutoff of zero. That is, we required a higher 
likelihood for a true miRNA than random background. This 
resulted in the prediction of 649 miRNAs in human (2,993 
total predictions, 33 1 known), 377 (3,063 total, 239 known) 
in chimpanzee, and 859 (3,538 total, 249 known) in rhesus 
macaque. Of these, 1 7 miRNAs were located in orthologous 
genomic regions in all three species and had a positive 
miRDeep score. One miRNA was predicted to be functional 
for both strands of the mature/star duplex in human. Seven 
miRNAs were independently described by others and were 
added to the new miRBase releases (version 1 7) during the 



Table 1 

Number of miRNAs Expressed in Each Tissue and in Each Species 
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preparation of this manuscript (Berezikov et al. 2006b; Goff 
et al. 2009; Creighton et al. 2010; Jima et al. 2010; 
Liao et al. 2010; Stark et al. 2010; Persson et al. 2011; 
Schotte et al. 2011; Dannemann et al. 2012). Twelve of 
the newly predicted miRNA candidates were tissue specific 
and eight of them were brain specific (supplementary 
fig. S2, Supplementary Material online). 

For the analyses presented in this paper, we used only the 
set of all miRBase-annotated miRNAs detected in our se- 
quencing data. For the between-species comparisons, we 
used miRNAs expressed in all three species and for the 
between-tissues comparisons, we used miRNAs expressed 
in all compared tissues in a given species. This is a rather strin- 
gent requirement, but it ensures that we compare bona fide 
miRNAs. 

miRNA Expression Differences 

We quantified miRNA expression and compared expression 
levels between species and tissues using the R package 
DESeq (Anders and Huber 2010). We found an approximately 
six times higher divergence in expression between tissues 
than between species. Divergence between tissues explained 
65% of the miRNA expression variance, whereas differences 
between species explained around 1 1 % (fig. 1c). Expression 
differences in brain explained the highest fraction of the tis- 
sue variance (41 %) followed by differences in heart (20%), 
liver (16%), kidney (14%,) and testis (9%) (fig. 1d). 

We compared the number of miRNAs with a significant 
expression difference between any two species (supplemen- 
tary fig. S3, Supplementary Material online). When compar- 
ing the fraction of miRNAs with significant difference in 
expression, we found that brain and heart show consistently 
fewer differences between species than other tissues. Be- 
tween humans and chimpanzees, 4% of all detected miRNAs 
differ significantly in expression in brain and 1 1 % in heart. In 
contrast, liver and testis had the most differentially 
expressed miRNAs (52% for both tissues). When comparing 
the fractions of differentially expressed miRNAs between 1) 
human and rhesus macaque and 2) chimpanzee and rhesus 
macaque, we found a significant difference between the 
two proportions for testis (45% and 25% for human- 
macaque and chimpanzee-macaque, respectively; P = 
0.006, Fisher's exact test) but no other tissue (P > 0.05). 

When we summed the number of differently expressed 
miRNAs over all three pairwise species comparisons, we ob- 
served a higher number of differently expressed miRNAs 
between human-macaque and between chimpanzee- 
macaque as compared with human-chimpanzee. The same 
pattern was consistently observed in all individual tissue com- 
parisons except in testis, which departed significantly from 
this trend. In testis, 38 miRNAs showed different expression 
between human and chimpanzee, 39 between chimpanzee 
and macaque, and 74 between human and macaque. 



Interestingly, miRNAs with expression differences be- 
tween tissues show often tissue-specific patterns, that is, 
these miRNAs differ in expression in only one tissue while 
expression is unchanged in the other four tissues. In hu- 
mans, 50% of these tissue-specific differentially expressed 
miRNAs were found in brain (chimpanzee 73%, rhesus ma- 
caque 43%). In total, we found that nine of the ten miRNAs 
with a consistent tissue-specific expression in all three 
species were brain specific. 

Correlation between miRNA Expression and Sequence 
Evolution 

To calculate the correlation between miRNA expression and 
sequence evolution, we used two measures: sequence con- 
servation among vertebrates and sequence variation within 
humans. We found positive Spearman correlations between 
average miRNA hairpin conservation score and miRNA ex- 
pression for all species and tissues (supplementary fig. S4, 
Supplementary Material online). All but one (chimpanzee tes- 
tis) were statistically significant (alpha = 0.05; supplementary 
table S3, Supplementary Material online). Conversely, we 
found negative Spearman correlations between the miRNA 
hairpin SNP number and miRNA expression though only in 
brain was the correlation statistically significant (supplemen- 
tary fig. S4 and table S3, Supplementary Material online). 

Functional Analysis of Differentially Expressed miRNAs 

To evaluate the effect of miRNA expression on the transcrip- 
tome, we used mRNA sequence tags quantified using se- 
quencing using the lllumina NG III Digital Gene Expression 
approach (Velculescu et al. 1995; Kircher et al. 2011). Of 
the 73 samples, 64 were common to both the mRNA and 
the miRNA studies (supplementary table S2, Supplementary 
Material online). mRNA targets of our expressed miRNAs 
were identified using TargetScanS (Lewis et al. 2003). 

We tested for functional enrichment of target genes in 
two of the GO (Ashburner et al. 2000) domains (molecular 
function and biological process) using FUNC (Prufer et al. 
2007). The multitissue/species comparisons allow us to dis- 
tinguish two groups of differentially expressed miRNAs. The 
first group is the tissue-specific miRNAs, which were iden- 
tified in the within-species comparisons (a total of 1 5 pair- 
wise tests: 3 species x 5 tissues). For a given species, we 
required that all miRNAs were expressed in all five tissues 
and filtered for miRNAs that were differentially expressed 
in one tissue compared with all others. The second group 
is the species-specific miRNAs, which were identified in 
the between-species comparisons (a total of 15 pairwise 
tests: 3 pairwise species comparisons x 5 tissues). For a given 
tissue, we required all miRNAs to be expressed in all species 
and filtered for those miRNAs that are differentially 
expressed in a species comparison. 
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We found 93 GO categories to be significantly enriched in 
at least one test for tissue-specific miRNA targets and 103 
categories in at least one test for species-specific miRNA tar- 
gets. The significantly enriched categories overlapped sub- 
stantially (62 categories) between the species- and the 
tissue-specific tests (supplementary table S4, Supplemen- 
tary Material online). We identified several highly connected 
clusters in the directed acyclic graph representation of the 
GO, which were related to development, cell communica- 
tion, gene and transcript expression, cell movement, regu- 
lation of metabolic, and biosynthetic processes and DNA/ 
protein/transcription binding (figs. 2a, 2b, and 3; supple- 
mentary figure S5, Supplementary Material online). 

We found that transcription factor-binding categories 
were significantly enriched — 23 of the total 30 GO tests 
performed. The signal was even more pronounced in the 
species-specific tests, where 1 3 of 1 5 tests were significant. 
Additionally, all tissue-specific tests in human were signifi- 
cant. In categories related to development, most of the 
significant enrichments were found in brain for both 
tissue-specific tests (16 of 25 significant tests in these cat- 
egories were brain specific) and species-specific compari- 
sons (15 of 21 significant tests in these categories were 
brain specific). Another cluster that showed enrichment 
for brain-related categories was cell communication in 
the species-specific test (of 22 significant tests, 13 were 
brain specific and 8 were testis specific). Additionally, we 
found that categories related to regulation of metabolic 
and biosynthetic processes (a total of 101 significant tests) 
were preferentially enriched in kidney (35) and liver (32) 
compared with brain (8), heart (6), and testis (10). 

Correlation between Transcription Factor Expression and 
miRNA Regulation 

Based on the functional enrichment for transcription factor 
activity, we sought to further investigate the regulatory 
relationship between miRNAs and transcription factors for 
the tissue comparisons. As a proxy for this functional 
relation, we used transcription factor expression and the 
number of miRNA-binding sites in each transcription factor. 
To obtain an annotation of transcriptions factors that is 
independent of GO categories, we designated a gene as 
a transcription factor if it was annotated in TRANSFAC 
(Kel et al. 2003). Using mRNA expression data and TRANS- 
FAC transcription factor annotation, we identified 94 tran- 
scription factors that were expressed in all five tissues. We 
identified differentially expressed transcription factors by 
performing pairwise tissue comparisons. We defined 
a transcription factor as tissue specific if it was differentially 
expressed between at least one tissue compared with the 
four others, that is, a transcription factor was allowed to 
show specific tissue expression patterns for one or more tis- 
sues. For each transcription factor, we obtained the number 



of 3' untranslated region (UTR)-binding sites for expressed 
miRNAs. We calculated the Spearman correlation (rho) be- 
tween the number of times a transcription factor is tissue- 
specific and the number of binding sites per expressed 
mRNA for an expressed miRNA and obtained a negative cor- 
relation (rho = -0.1 5, P = 0. 1 6). We compared the distri- 
butions of the number of binding sites for expressed miRNAs 
between tissue-specific and nontissue-specific transcription 
factors and found a nonstatistically significant shift between 
the distributions (supplementary fig. S6, Supplementary 
Material online). We repeated the two analyses using only 
binding sites for miRNAs that are highly differentially ex- 
pressed between tissues (tissue specific for more than 
two tissues) and obtained a negative correlation (rho = 
-0.18, P = 0.09) and a significant difference between 
the number of miRNA-binding sites for the two groups of 
transcription factors (Wilcoxon rank sum test, P = 0.03). 
Transcription factors with fewer binding sites for expressed 
miRNAs tend to show a higher variability in their expression 
between tissues than those with more binding sites. 

Correlation of Expression between miRNAs and Their 
Predicted Target Genes 

To determine the functional relationship between miRNAs 
and their predicted target genes, we correlated the expres- 
sion values of miRNAs and genes. For each miRNA-mRNA 
target pair, we required that both the miRNA and the mRNA 
were expressed in all species or tissues. We then carried out 
two tests to determine the level of correlation at different 
levels of granularity. Each test calculates the Pearson corre- 
lation coefficient (r) for each pair. 

For our first test, we correlated miRNA and mRNA expres- 
sion between tissues in each species. We averaged the 
miRNA and mRNA expression values over all individuals 
per tissue and species. We then calculated the Pearson cor- 
relation coefficient between miRNA expression and mRNA 
expression levels for the five tissues. The fraction of negative 
correlations observed was not different from the fraction 
obtained by randomly assigning miRNAs and mRNAs. When 
comparing genes with one miRNA-binding site to genes 
with more than one, we found a significant signal in all spe- 
cies (to obtain significance, we randomly permuted the 
number of binding sites in genes, p_hu < 0.01, p_ch < 
0.01, p_rh < 0.01; fig. 4). However, when we restricted 
our analysis to genes with only one miRNA-binding site, 
we found that the fraction of negative correlations was sig- 
nificantly higher than by randomly assigning miRNAs and 
mRNAs in human but none of the other species (p_hu = 
0.02, p_ch = 0.26, p_rh = 0.65; fig. 4). 

For our second test, we sought to test the correlation of 
miRNA to target mRNA between species. For each tissue and 
miRNA-mRNA pair, we arrived at three datapoints that are 
used as input for the individual correlations. The correlations 
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were calculated for each pair, and the resulting fraction of 
negative correlations was tested against random target 
assignments. We also tested the fraction of negative correla- 
tions dependent on the number of miRNA-binding sites in the 



genes. None of the tissues gave a significant excess of neg- 
ative correlations, and the amount of negative correlation 
was not different from the amount obtained from random 
assignments of genes to miRNAs (all P values > 0.05). The 
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excess of negative correlation also did not differ depending 
on the number of miRNA-binding sites in the genes. 

It has been proposed that miRNAs exert their effect by 
setting the mean and reducing the variance of the genes 



that they regulate, and in this way, they stabilize phenotypes 
a process called canalization (Wu et al. 2009). Our results 
are consistent with miRNAs setting the mean expression 
of the genes they regulate. For technical reasons, we are 
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not able to test where there is also a reduction in the 
variance of mRNA expression. We found a correlation 
between the mean expression level of the genes and the 
relative variance to the mean in the mRNA data set, that 
is, miRNAs tend to target highly expressed genes, and these 
genes also show a higher variance due to the fact that their 
expression level and their variance are not independent. 
Unfortunately, we found no general normalization that 



eliminated the correlation between the variance estimate 
and the expression of a given gene. 

Discussion 

In this study, we have characterized small RNA libraries in 
five tissues of humans, chimpanzees, and rhesus macaques 
using high-throughput sequencing. For four of the five 
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tissues, the majority of the sequencing reads mapped to 
known miRNAs. In contrast, testis showed only a small per- 
centage (9%) of reads matching known miRNAs. A larger 
fraction mapped to piRNAs (17%) while about one-third 
mapped to the genome but did not overlap any known 
RNA annotation. PiRNAs are known to be expressed specif- 
ically in germ line and gonadal somatic cells (Siomi et al. 
2011) and are associated with the PIWI proteins, which 
are indispensable proteins for germ line development (Siomi 
et al. 201 1) in many animals (Ghildiyal and Zamore 2009). 
They are divided in two classes: the first class is involved in 
silencing transposons, whereas the function of the second 
class remains unknown (Aravin et al. 2007). PiRNAs have 
a length distribution of between 25 and 35 nt (Aravin et al. 
2001). Our library preparation protocol aims at extracting 
molecules of up to 25 nt in length, implying that the set 
of piRNAs discovered in our study is biased toward smaller 
sizes. Despite this limitation, we were able to detect a bias 
for 5' uridine (U) residue (Malone and Hannon 2009) and an 
excess of adenosine at position 10 (Friedlander et al. 2009) 
that characterize piRNAs. Interestingly, testis reads that 
mapped to each species' corresponding genome but did 
not overlap any known RNA annotation showed the same 
patterns, suggesting that a large fraction of piRNAs remain 
to be characterized, as they are not among the —32,000 
human piRNAs reported in the database we used for map- 
ping (Pang et al. 2007). 

Summing over all tissues, we were able to detect approx- 
imately 80% of all known miRNAs for each of the three spe- 
cies. When counting the number of miRNAs expressed, we 
observed a difference between tissues. In all three species, 
brain and testis show the most diverse miRNA repertoire. For 
testis, this result is surprising given the reduced power that 
we have to detect lowly expressed miRNAs because of the 
large fraction of reads representing other small RNAs that 
may or may not be transcribed from unannotated parts 
of the genome. In contrast, the miRNA repertoires of heart 
and liver are less diverse (table 1). Given that miRNAs have 
been shown to regulate mRNA abundance it is to be ex- 
pected that the diversity of expressed miRNAs is linked to 
the number of different mRNAs expressed in each tissue. 
MRNA expression diversity was shown in previous studies 
to be highest in testis (Khaitovich et al. 2005; Kircher 
et al. 2011). However, gene expression diversity in brain 
was similar to the diversity seen in liver, kidney, and heart, 
and not more diverse, as observed for miRNA expression in 
the brain. This discrepancy between mRNA and miRNA di- 
versity is further supported by our observation that a large 
fraction of the newly predicted miRNAs were brain specific 
(supplementary fig. S2, Supplementary Material online). Our 
data thus suggest that a larger fraction of genes may be reg- 
ulated by miRNAs in brain as compared with other tissues. If 
true, this would lend further support to the central role at- 
tributed to regulatory RNAs in brain function (Mattick 



201 1). However, we also note that brain is an intensively 
studied tissue and that miRNAs annotated in miRBase could 
be potentially biased toward brain-specific miRNAs. A bias in 
annotation cannot, however, explain why newly discovered 
miRNAs are often brain specific. MiRNAs show high se- 
quence conservation between species (Pasquinelli et al. 
2000; Hertel et al. 2006; Chen and Rajewsky 2007; Stark 
et al. 2007). In this study, we found a similar high conser- 
vation in their expression between primate species. In agree- 
ment with previous gene expression studies, the divergence 
between tissues is higher than the divergence between spe- 
cies. This is reflected in the high amount of variance (65%) 
explained by the factor "tissue" in the linear model that we 
applied (fig. 1 c). The high conservation of miRNA expression 
levels and tissue-specific expression patterns were previously 
observed in a comparative study of 26 tissues in humans and 
rodents (Landgraf et al. 2007). In our study, miRNAs ex- 
pressed in brain explained the highest fraction of tissue var- 
iance (fig. 1 d). It has been shown that brain has the smallest 
mRNA expression differences between species, and it was 
concluded that purifying selection is extremely efficient at 
eliminating mutations that modify gene expression in this 
tissue (Khaitovich et al. 2005, 2006). We found that many 
miRNAs are specifically expressed in brain. We therefore hy- 
pothesize that miRNAs, in addition to purifying selection, 
may explain the highly conserved mRNA expression patterns 
for this tissue. 

Although the expression differences in miRNAs between 
species are small, they are sufficient to pinpoint miRNAs 
with a difference in expression between the primate species 
studied. Humans and chimpanzees show more similarity 
than either shows to rhesus macaques. The differences be- 
tween species are particularly small in brain and heart and 
much bigger in liver and testis. It has been previously shown 
that mRNA expression is subjected to different levels of con- 
straint in different tissues (Khaitovich et al. 2005, 2006). Us- 
ing rhesus macaque as outgroup to assign differences, we 
found no excess of differential expression on the human or 
chimpanzee lineage for the tissues brain, kidney, liver, and 
heart. In testis, however, we observed that chimpanzees and 
rhesus macaques are much closer to one another than either 
is to humans (supplementary fig. S3, Supplementary Mate- 
rial online). Reproductive pressures such as sperm competi- 
tion may be similar for chimpanzees and rhesus macaques 
(Dixson 1998), and it has been hypothesized that selection 
may have shaped both the protein sequence evolution and 
the mRNA expression of male reproductive genes (Wyckoff 
et al. 2002; Khaitovich et al. 2005, 2006). Our data suggest 
that adaptive changes may have shaped the expression pat- 
terns of miRNAs in testis on the human lineage. 

We sought to study the relation between expression and 
sequence evolution. We found a positive correlation be- 
tween miRNA sequence conservation and expression level. 
It has previously been shown in humans and Drosophila that 
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purifying selection is weaker in lowly expressed miRNAs and 
that they therefore have less constraint in sequence evolu- 
tion than highly expressed miRNAs (Lu et al. 2008; Liang and 
Li 2009). Highly expressed miRNAs are important regulators 
of gene expression, and their sequence conservation is 
therefore accordingly higher. The levels of purifying selec- 
tion can be also assessed by the within species diversity. 
We took advantage of the extensive human diversity data- 
bases and found negative correlations between miRNA 
expression and the number of miRNA polymorphisms. 
A similar effect has been reported in humans (Liang and 
Li 2009). However, the effect was not as strong as the cor- 
relation between expression and conservation. The lack of 
polymorphism data available for other primate species 
meant that we were not able to perform the same analysis 
on these. 

In addition to known miRNAs, we predicted 17 new 
miRNAs in the three species studied. It has been debated 
whether newly predicted miRNAs are functional or not, es- 
pecially because they tend to have lower expression and 
more sequence divergence than known miRNAs (Lu et al. 
2008; Liang and Li 2009). While this paper was under re- 
view, 7 of our 17 newly predicted miRNAs were reported 
by other studies and add to miRBase (supplementary 
fig. S2, Supplementary Material online). The fact that these 
miRNAs were independently predicted by other groups, to- 
gether with their orthologous location in the three primate 
species and their tissue-specific expression patterns support 
that our predictions are reliable (supplementary table S5, 
Supplementary Material online). Our comparison of the 
expression levels of newly predicted miRNAs with known 
miRNAs showed that the former group tended to have 
lower expression values (supplementary fig. S7, Supplemen- 
tary Material online). This is in agreement with previous 
studies and supports the hypothesis that newly emerged 
miRNAs are raw material for evolution that only have 
weak or negligible impact on gene expression after their 
emergence (Liang and Li 2009). 

Mammalian miRNAs exert their regulatory action by de- 
creasing target mRNA levels (Guo et al. 2010). We there- 
fore expected to find a negative correlation between 
miRNAs and their target genes. We used miRNA targets 
that were computed based on genome-wide prediction al- 
gorithms that take sequence complementarity into ac- 
count. For these predicted targets, we performed the 
analysis at tissue and species level. We did not find an ex- 
cess of negative correlation between miRNA and mRNA ex- 
pression in the between tissues and between species 
comparisons. However, when we restricted our set of tar- 
gets to genes with only one binding site for expressed 
miRNAs, we did find an excess of negative correlation in 
the between tissues comparison, whereas the between 
species comparison gave no signal. This difference in 
results may be due to a difference in power; tissue diver- 



gence exceeds species divergence, with tissue differences 
explaining the majority of variance in miRNA expression 
(fig. 1c). It is possible that technical variation exceeds 
the effect exerted by miRNAs at the gene level in the 
between species comparisons, whereas the comparison 
between tissues yields significant results due to more pro- 
nounced differences in expression. The lack of power in the 
between-species comparison may be further exacerbated 
by the use of human-based target prediction databases 
that could potentially obscure species-specific effects. 
However, a substantial source of error may lie in the high 
false positive rate of computational miRNA target predic- 
tion, which poses a challenge for understanding the impact 
of miRNAs on gene expression regulation (Alexiou et al. 
2009). Another source of noise for finding the regulatory 
relation between the expression of miRNA-mRNA pairs is 
the involvement of other regulatory molecules that hinder 
a clear measurement of the effect of miRNAs alone. Tran- 
scription factors, in particular, can regulate transcription 
positively and negatively (Hobert 2008) and their effect 
strength has been reported to be larger than that of 
miRNAs: in Caenorhabditis elegans individual deletion of 
miRNA loci only caused developmental and morphological 
defects in <10% of the cases (Miska et al. 2007), whereas 
RNA interference experiments for transcription factors re- 
sulted in observable effects in —30% of the cases reviewed 
in Hobert (2008). The action of transcription factors may 
therefore lead to false signals of apparent positive correla- 
tion between miRNA-mRNA pairs, and mRNA expression 
differences have to be seen as the result of the sum of 
many regulating factors impacting the expression level. 

The pronounced excess of negative correlations, we ob- 
serve when using only genes with fewer binding sites for 
expressed miRNAs indicates that the power to measure 
the functional relation between miRNA-mRNA pairs is 
higher for this group of genes. Genes regulated by multiple 
expressed miRNAs introduce more noise and hinder our 
ability to measure the impact of each miRNA on mRNA 
expression values. Two studies of development and aging 
in humans and macaque brains (Somel et al. 2010), and 
humans, chimpanzee, and macaque brains (Somel et al. 
201 1) have shown that age-related miRNA expression pro- 
files are more negatively correlated with their targets than 
random expectations. A similar effect has been reported by 
comparing mRNAs, miRNAs, and proteins from human and 
chimpanzees using two different brain regions (Hu et al. 
2011). 

By comparing multiple species and tissues, we found that 
genes targeted by differentially expressed miRNAs were en- 
riched in some GO categories, whereas genes targeted by 
miRNAs with uniform expression showed no signal of en- 
richment. The significant categories formed functionally re- 
lated clusters in the GO graphs for the biological process and 
molecular function domains (fig. 2a and b). A significant 
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enrichment for genes with transcription factor-binding ac- 
tivity was found consistently in the majority of tests both 
between species and between tissues (fig. 3). This suggests 
that miRNAs, as downregulators of gene expression, prefer- 
entially modulate transcription factors, which are them- 
selves transacting regulatory molecules. It has been 
shown that transcription factors that are differentially ex- 
pressed between human and chimpanzee brains preferen- 
tially regulate genes involved metabolism and transcription 
among others (Nowick et al. 2009). Additionally, in plants, 
miRNAs preferentially regulate transcription factors that are 
involved in development (Rhoades et al. 2002; Jones- 
Rhoades and Bartel 2004). We speculate that miRNAs 
may therefore also have an indirect influence on the expres- 
sion of these types of genes. 

As a second line of evidence for the connection between 
miRNAs and transcription factors, we measured the func- 
tional relationship between transcription factor expression 
differences and miRNA regulation between tissues. To do that 
we used as a proxy the number of binding sites for expressed 
miRNA present in transcription factors 3' UTRs. We found 
that differentially expressed transcription factors contain 
fewer miRNA-binding sites than transcription factors that 
are not differentially expressed. MiRNAs have been specu- 
lated to have the function of stabilizing the expression of their 
targets (Kircher et al. 2008; Wu et al. 2009). In our between- 
tissue comparison, this would mean that the miRNA expres- 
sion level is adjusted to maintain homogeneous mRNA 
expression. Transcription factors that are differentially ex- 
pressed show a depletion in miRNA-binding sites. We hypoth- 
esize that these transcription factors achieve differential 
expression by reducing the number of miRNA-binding sites 
and therefore avoiding the control by expressed miRNAs. 

Further studies integrating miRNA-mRNA profiles with 
data sets produced by shotgun proteomics and ribosome 
profiling will further improve the understanding of gene ex- 
pression and gene expression regulation in primates and 
help to disentangle the relative contributions of different 
gene expression regulatory machinery and their impact 
on phenotype in primate evolution. 

Supplementary Material 

Supplementary figures S1-S7 and tables S1-S5 are available 
at Genome Biology and Evolution online (http:// 
www.gbe.oxfordjournals.org/). 

Acknowledgments 

We thank Thomas Giger for performing the frozen tissue 
dissections; Ines Drinnenberg, Matthias Meyer, and the Se- 
quencing Group of the Max Planck Institute for Evolutionary 
Anthropology for providing sequencing runs; Martin Kircher 
for technical assistance with sequence processing; Marike 
Schreiber for assistance with the figure preparation; Cesare 



de Filippo, Martin Kircher, Michael Lachmann, Mehmet 
Somel, and Martin Surbeck for useful discussions during 
the project and for comments on the manuscript. The pro- 
ject was funded by a grant of the Max Planck Society. 

Literature Cited 

Alexiou P, Maragkakis M, Papadopoulos GL, Reczko M, 
Hatzigeorgiou AG. 2009. Lost in translation: an assessment and 
perspective for computational microRNA target identification. 
Bioinformatics 25:3049-3055. 

Anders S, Huber W. 2010. Differential expression analysis for sequence 
count data. Genome Biol. 1 1 :R1 06. 

Aravin AA, Hannon GJ, Brennecke J. 2007. The Piwi-piRNA pathway 
provides an adaptive defense in the transposon arms race. Science 
318:761-764. 

Aravin AA, et al. 2001. Double-stranded RNA-mediated silencing 
of genomic tandem repeats and transposable elements in the 
D. melanogaster germline. Curr Biol. 1 1:1017-1027. 

Ashburner M, et al. 2000. Gene ontology: tool for the unification of 
biology. The Gene Ontology Consortium. Nat Genet. 25:25-29. 

Baek D, et al. 2008. The impact of microRNAs on protein output. Nature 
455:64-71. 

Baev V, Daskalova E, Minkov I. 2009. Computational identification of 
novel microRNA homologs in the chimpanzee genome. Comput Biol 
Chem. 33:62-70. 

Berezikov E, et al. 2005. Phylogenetic shadowing and computational 

identification of human microRNA genes. Cell 120:21-24. 
Berezikov E, et al. 2006a. Diversity of microRNAs in human and 

chimpanzee brain. Nat Genet. 38:1375-1377. 
Berezikov E, et al. 2006b. Many novel mammalian microRNA candidates 

identified by extensive cloning and RAKE analysis. Genome Res. 

16:1289-1298. 

Brameier M. 2010. Genome-wide comparative analysis of microRNAs in 

three non-human primates. BMC Res Notes. 3:64. 
Chen K, Rajewsky N. 2007. The evolution of gene regulation by 

transcription factors and microRNAs. Nat Rev Genet. 8:93-103. 
Creighton CJ, Reid JG, Gunaratne PH. 2009. Expression profiling of 

microRNAs by deep sequencing. Brief Bioinformatics 10:490-497. 
Creighton CJ, et al. 2010. Discovery of novel microRNAs in female 

reproductive tract using next generation sequencing. PLoS One 5:e9637. 
Dannemann M, Nickel B, Lizano E, Burbano HA, Kelso J. 2012. 

Annotation of primate miRNAs by high-throughput sequencing of 

small RNA libraries. BMC Genomics 13:116. Advance Access 

publication March 27, 2012. doi:1 0.1 186/1471-21 64-13-1 1 6. 
Dixson AF. 1998. Primate sexuality. Oxford: Oxford University Press. 
Erson AE, Petty EM. 2008. MicroRNAs in development and disease. Clin 

Genet. 74:296-306. 
Friedlander MR, et al. 2008. Discovering microRNAs from deep 

sequencing data using miRDeep. Nat Biotechnol. 26:407-415. 
Friedlander MR, et al. 2009. High-resolution profiling and discovery of 

planarian small RNAs. Proc Natl Acad Sci USA. 106:11546-11551. 
Friedman RC, Farh KK, Burge CB, Bartel DP. 2009. Most mammalian 

mRNAs are conserved targets of microRNAs. Genome Res. 19:92-105. 
Ghildiyal M, Zamore PD. 2009. Small silencing RNAs: an expanding 

universe. Nat Rev Genet. 10:94-108. 
Girard A, Sachidanandam R, Hannon GJ, Carmell MA. 2006. 

A germline-specific class of small RNAs binds mammalian Piwi 

proteins. Nature 442:199-202. 



Genome Biol. Evol. 4(4):552-564. doi:10.1093/gbe/evs033 Advance Access publication March 27, 2012 



563 



Dannemann et al. 



GBE 



Goff LA, et al. 2009. Ago2 immunoprecipitation identifies predicted 

microRNAs in human embryonic stem cells and neural precursors. 

PLoS One 4:e7192. 
Griffiths-Jones S, Saini HK, van Dongen S, Enright Al. 2008. miRBase: 

tools for microRNA genomics. Nucleic Acids Res. 36:D1 54-D1 58. 
Guo H, Ingolia NT, Weissman JS, Bartel DP. 2010. Mammalian 

microRNAs predominantly act to decrease target mRNA levels. 

Nature 466:835-840. 
Hertel J, et al. 2006. The expansion of the metazoan microRNA 

repertoire. BMC Genomics 7:25. 
Hinrichs AS, et al. 2006. The UCSC Genome Browser Database: update 

2006. Nucleic Acids Res. 34:D590-598. 
Hobert O. 2008. Gene regulation by transcription factors and micro- 
RNAs. Science 319:1785-1786. 
Hu HY, et al. 2011. MicroRNA expression and regulation in human, 

chimpanzee, and macaque brains. PLoS Genet. 7:e1 002327. 
Jima DD, et al. 2010. Deep sequencing of the small RNA transcriptome 

of normal and malignant human B cells identifies hundreds of novel 

microRNAs. Blood 1 1 6:e1 1 8-e1 27. 
Jones-Rhoades MW, Bartel DP. 2004. Computational identification of 

plant microRNAs and their targets, including a stress-induced 

miRNA. Mol Cell. 14:787-799. 
Kawaji H, et al. 2008. Hidden layers of human small RNAs. BMC 

Genomics 9:1 57. 

Kel AE, et al. 2003. MATCH: a tool for searching transcription factor 

binding sites in DNA sequences. Nucleic Acids Res. 31:3576-3579. 
Khaitovich P, et al. 2005. Parallel patterns of evolution in the genomes and 

transcriptomes of humans and chimpanzees. Science 309: 1 850-1 854. 
Khaitovich P, et al. 2006. Positive selection on gene expression in the 

human brain. Curr Biol. 1 6:R356-R358. 
King MC, Wilson AC. 1975. Evolution at two levels in humans and 

chimpanzees. Science 188:107-116. 
Kircher M. 2011. Understanding and improving high-throughput 

sequencing data production and analysis [dissertation]. [Leipzig 

(Germany)]: Fakultat fur Mathematik und Informatik. 
Kircher M, Bock C, Paulsen M. 2008. Structural conservation versus 

functional divergence of maternally expressed microRNAs in the 

Dlk1/Gtl2 imprinting region. BMC Genomics 9:346. 
Kircher M, Stenzel U, Kelso J. 2009. Improved base calling for the 

lllumina Genome Analyzer using machine learning strategies. 

Genome Biol. 10:R83. 
Landgraf P, et al. 2007. A mammalian microRNA expression atlas based 

on small RNA library sequencing. Cell 129:1401-1414. 
Lee RC, Feinbaum RL, Ambros V. 1993. The C. elegans heterochronic 

gene lin-4 encodes small RNAs with antisense complementarity to 

lin-14. Cell 75:843-854. 
Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB. 2003. 

Prediction of mammalian microRNA targets. Cell 115:787-798. 
Liang H, Li WH. 2009. Lowly human expressed microRNA genes evolve 

rapidly. Mol Biol Evol. 26:1195-1198. 
Liao JY, et al. 2010. Deep sequencing of human nuclear and cytoplasmic 

small RNAs reveals an unexpectedly complex subcellular distribution 

of miRNAs and tRNA 3' trailers. PLoS One 5:e10563. 
Lim LP, et al. 2005. Microarray analysis shows that some microRNAs 

downregulate large numbers of target mRNAs. Nature 433:769-773. 
Lu J, et al. 2008. Adaptive evolution of newly emerged micro-RNA 

genes in Drosophila. Mol Biol Evol. 25:929-938. 
Malone CD, Hannon GJ. 2009. Small RNAs as guardians of the genome. 

Cell 136:656-668. 
Mattick JS. 2011. The central role of RNA in human development and 

cognition. FEBS Lett. 585:1600-1616. 



Miska EA, et al. 2007. Most Caenorhabditis elegans microRNAs are 
individually not essential for development or viability. PLoS Genet. 
3:e215. 

Nowick K, Gernat T, Almaas E, Stubbs L. 2009. Differences in human 
and chimpanzee gene expression patterns define an evolving 
network of transcription factors in brain. Proc Natl Acad Sci USA. 
106:22358-22363. 

Pang KC, et al. 2007. RNAdb 2.0— an expanded database of 
mammalian non-coding RNAs. Nucleic Acids Res. 35:D178-D182. 

Pasquinelli AE, et al. 2000. Conservation of the sequence and temporal 
expression of let-7 heterochronic regulatory RNA. Nature 408: 
86-89. 

Persson H, et al. 201 1 . Identification of new microRNAs in paired normal 

and tumor breast tissue suggests a dual role for the ERBB2/Her2 

gene. Cancer Res. 71:78-86. 
Prufer K, et al. 2007. FUNC: a package for detecting significant 

associations between gene sets and ontological annotations. BMC 

Bioinformatics 8:41 . 
Prufer K, et al. 2008. PatMaN: rapid alignment of short sequences to 

large databases. Bioinformatics 24:1530-1531. 
Rhoades MW, et al. 2002. Prediction of plant microRNA targets. Cell 

110:513-520. 

Schotte D, et al. 201 1. Discovery of new microRNAs by small RNAome 

deep sequencing in childhood acute lymphoblastic leukemia. 

Leukemia 25:1389-1399. 
Selbach M, et al. 2008. Widespread changes in protein synthesis 

induced by microRNAs. Nature 455:58-63. 
Sherry ST, et al. 2001. dbSNP: the NCBI database of genetic variation. 

Nucleic Acids Res. 29:308-311. 
Siepel A, et al. 2005. Evolutionarily conserved elements in vertebrate, 

insect, worm, and yeast genomes. Genome Res. 15:1034-1050. 
Siomi MC, Sato K, Pezic D, Aravin AA. 2011. PlWI-interacting small 

RNAs: the vanguard of genome defence. Nat Rev Mol Cell Biol. 

12:246-258. 

Somel M, et al. 2010. MicroRNA, mRNA, and protein expression link 
development and aging in human and macaque brain. Genome Res. 
20:1207-1218. 

Somel M, et al. 2011. MicroRNA-driven developmental remodeling in 
the brain distinguishes humans from other primates. PLoS Biol. 
9:e1001214. 

Stark A, et al. 2007. Discovery of functional elements in 12 
Drosophila genomes using evolutionary signatures. Nature 
450:219-232. 

Stark MS, et al. 2010. Characterization of the melanoma miRNAome by 

deep sequencing. PLoS One 5:e9685. 
Velculescu VE, Zhang L, Vogelstein B, Kinzler KW. 1995. Serial analysis 

of gene expression. Science 270:484-487. 
Wightman B, Ha I, Ruvkun G. 1993. Posttranscriptional regulation of the 

heterochronic gene lin-14 by lin-4 mediates temporal pattern 

formation in C. elegans. Cell 75:855-862. 
Wu CI, Shen Y, Tang T. 2009. Evolution under canalization and the dual 

roles of microRNAs: a hypothesis. Genome Res. 19:734-743. 
Wyckoff GJ, Li J, Wu CI. 2002. Molecular evolution of functional 

genes on the mammalian Y chromosome. Mol Biol Evol. 

19:1633-1636. 

Yue J, Sheng Y, Orwig KE. 2008. Identification of novel homologous 
microRNA genes in the rhesus macaque genome. BMC Genomics 9:8. 



Associate editor: George Zhang 



564 Genome Biol. Evol. 4(4):552-564. doi:10.1093/gbe/evs033 Advance Access publication March 27, 2012 



